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A crystal of emerging magnetic charges is expected in the phase diagram of the dipolar kagome spin ice. An 
observation of charge crystallites in thermally demagnetized artificial spin ice arrays has been recently 
reported by S. Zhang and coworkers 1 and explained through the thermodynamics of the system as it 
approaches a charge- ordered state. Following a similar approach, we have generated a partial order of 
magnetic charges in an artificial kagome spin ice lattice made out of ferrimagnetic material having a Curie 
temperature of 475 K. A statistical study of the size of the charge domains reveals an unconventional 
sawtooth distribution. This distribution is in disagreement with the predictions of the thermodynamic 
model and is shown to be a signature of the kinetic process governing the remagnetization. 

Frustrated spin ices are usually described by Hamiltonians containing only spin-spin interactions. These 
interactions along with the lattice topology lead to the so called spin ice rules (also named Bernal-Fowler 
rules). In a 2D kagome lattice composed of vertices linking 3 spins (figure la), the spin ice rule takes the 
following form: allowed states are either made of 1 incoming and 2 outgoing spins (top raw of figure lb) or by their 
equivalents obtained by a time reversal transformation (i.e. 1 outgoing and 2 incoming spins, bottom raw of 
figure lb). The spin balance at a vertex i can naturally be represented as a magnetic charge Q f . This entails a 
simplified spin ice rule: only ± 1 charge at each vertex are allowed in the ice state (figure lb). Therefore, a violation 
of this rule implies the presence of a +3 or —3 charge (figure lc). 

By explicitly including long-range spin-spin interactions, the kagome dipolar spin ice model (DSI) exhibits a 
richer behavior than the nearest neighbor interaction version, while it still preserves the ice rule 2 " 4 . Indeed, the 
variation of the DSI entropy with respect to the system temperature shows four plateaux each of them being 
associated to a different phase (figure lg). At high temperatures, the DSI is in a paramagnetic regime violating the 
ice rule. When reducing the temperature, the system enters into the spin ice manifold, where the ice rule is 
unanimously obeyed by all vertices (i.e. composed only of allowed states). However, in this regime spins and 
magnetic charges are not showing an explicit order. Going to lower temperatures, the system reaches a state 
named spin ice II. In this recently predicted regime 3 ' 4 , while spins are still allowed to fluctuate, the magnetic charge 
distribution is frozen. And, even though magnetic charges do not explicitly appear in the Hamiltonian, a crystal of 
emerging charges is formed (figure Id and le). Parts of this crystal, made of a succession of positive and negative 
emerging charges, have been recently experimentally observed by Zhang and col. 1 . Finally, at low temperatures, 
the system reaches its ground state which is characterized by a spin-loop crystal in which both the spins and the 
magnetic charges are ordered (figure If). 

However, experimentally reaching the spin ice II manifold in artificial systems remains challenging, since the 
effective temperature has to be controlled and strongly reduced. Different field procedures have been used to 
explore the energy manifold of artificial spin ice systems (particularly square ice) 5 " 7 . Concerning the Kagome 
lattice, it has been shown that the DSI model has to be considered 8 but no spin ice II have been observed 8 ' 9 . Purely 
thermal demagnetization was only recently reported 110 " 13 , without quantitative comparison to DSI for Kagome 
ice. 

In this study, we use a Co 0 . 7 Gdo.3 ferrimagnetic alloy thin film to form the kagome lattice. This film has a 
magnetization of 350 kA/m and a Curie temperature (T c ) of 475 K (see supplementary material A). Heating the 
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Figure 1 | (a) AFM image of CoGd connected nanomagnets in the kagome geometry, (b) Sketches of the six spin configurations allowed by the ice rule at a 
vertex. The blue and red spheres of small radius correspond respectively to the -1 and + 1 magnetic charges located at the vertex, (c) Sketches of the two 
forbidden states. The blue and red spheres of large radius correspond to the —3 and +3 magnetic charges, (d) Spin and charge configurations in the Spin 
Ice II state, (e) A charge crystal in the Spin Ice II state, (f) A spin loop crystal in the ground state, (g) Entropy and charge correlator as a function of 
temperature in the kagome DSI framework. The temperature is normalized to the nearest neighbor coupling coefficient. 



film up to 500 K does influence neither its magnetic properties nor 
its amorphous structure (as checked by X-ray diffraction). This 
material thus provides the possibility of a thermal demagnetization 
by cooling the sample from its Curie point without inflicting any 



damage to the system. The studied arrays are formed of 342 con- 
nected nanomagnets, arranged in a kagome geometry (see figure la). 
The distance between adjacent vertices is 500 nm. At room temper- 
ature, due to the shape anisotropy, each element presents a uniform 
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Figure 2 | (a) XMCD-PEEM image, (b) Spin configuration deduced from the XMCD-PEEM image and (c) associated magnetic charge 
configuration along with domains of emerging magnetic charges (green and yellow) after a thermal demagnetization, (d) Distribution of the measured 
charge correlator (in red) and calculated distribution in the SI model (black curve), (e) Zoom on the array to show spins, emerging charges and 
charge domains. 



magnetic state, with magnetization oriented along the element's long 
axis. This mimics a multi-axis Ising spin system. 

In a first step, the array is magnetized by the application of an in- 
plane magnetic field. After reducing the external field to zero, the 
array is brought above its Curie point and therefore beyond the 
magnetic order temperature of the alloy. Subsequently, the sample 
is cooled down. The magnetic configuration of our kagome struc- 
tures was experimentally determined by XMCD-PEEM, a direct 
imaging technique allowing us to image the orientation of the mag- 
netization vector in each nano -magnet. A sample image is presented 
in figure 2a, along with the corresponding scheme of the spin con- 
figuration (figure 2b). After cooling below T a every magnet of the 
network is again uniformly magnetized, but the memory of the satu- 
rated state is lost. Analysis of the "spins" and vertices distribution 
shows that the demagnetization is effective with a remanent magnet- 
ization of 3.2% (see supplementary materials C). As the spin ice II 
state is associated with a charge order, the emerging magnetic 
charges at each vertex were calculated (figure 2c-e). 

Our data indicate that the ice rule is strictly obeyed among more 
than 3800 observed vertices. Consequently, only ±1 charges are 
associated to these vertices. In order to get an insight about the charge 
order, the nearest-neighbor charge correlator <Q,-.Qi+i> has been 
computed for 18 different kagome networks. As shown in figure 2d, 
the charge correlator distribution differs clearly from the one 
expected for the spin ice state 8 . We found an average value of 
—0.22 lower than —1/9, the <Q/.Q/+i> value in the spin ice state 
(c.f. figure lg). Nevertheless, it remains higher than —1, the value 
characterizing a fully ordered spin ice II state (perfect crystal of 
emerging charges). However, quantitative analysis of the spatial 
charge distribution reveals the presence of a local charge order. 

Figure 2c presents ordered charge areas that belong to antiphase 
charge domains, highlighted in yellow and green. The size of these 
domains can extend up to 80 vertices i.e. 37% of the total number of 
vertices in a network. The domains are delimited by domain walls 
separating vertices of equal charges (figure 2e). These domain walls 



are analogous to antiphase boundaries. The spin configurations 
within the charge domains do not exhibit trivial spin order (neither 
saturated nor crystal loop state). 

The domain size distribution extracted from the 18 studied net- 
works is presented in figure 3a. The general trend is a decrease of the 
number of domains with their size (inset in figure 3a). However, a 
closer look at the small domain sizes reveals a non-monotonic dis- 
tribution: domains of 1, 3 and 5 vertices appear to be underrepre- 
sented. This sawtooth shape in the domain size distribution is 
unexpected. We can exclude here any artifact due to our sample as 
a "classic" rotational field- decaying demagnetization leads to a 
monotonous domain size distribution very similar to the one 
expected for the spin ice phase (figure 3b). 

Discussion 

To establish the origins of this astonishing domain size distribution, 
Monte Carlo simulations were employed to study the transition 
towards the observed partial order of charges admitting thermodyn- 
amic equilibrium conditions are met. Figure 3c represents the evolu- 
tion of the number of domains (with size from 1 to 8 vertices) as a 
function of the charge correlator (cf. methods). For these domain 
sizes, the number of domains monotonically decreases with the 
charge correlator (and temperature) in favor of bigger domains. 
For charge correlators above —0.55, the number of domains per 
array monotonically decreases with the size of the domain. Below 
— 0.55, the number of 2 vertices charge domains exceeds the one of 
size 1 (red/black dot on figure 3c). The number of domains contain- 
ing three vertices always remains higher than the number of domains 
having 4 or more vertices. Also note that the occurrence of size 5 and 
6 domains (as well as size 7 compared to 8), are quite similar for a 
nearest neighbor charge correlator < — 0.55 (The difference is below 
0.1 domain/network). This behavior thus contrasts with our experi- 
mental findings. In figure 3a, the experimentally observed domain 
size distribution is directly compared to the one expected for the DSI 
model at equilibrium for the same charge correlator. The number of 
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Figure 3 | (a) Distribution of domain sizes (in vertices) after a thermal demagnetization and comparison to the DSI model (error bars represent the 
fluctuation of the model for the number of samples considered in our experiments), (b) Distribution of domain sizes after a rotating field 
demagnetization and comparison to the spin ice model, (c) and (d) Numbers of domains per array for domains including 1 to 8 vertices in the DSI 
framework and in the thermo-kinetic framework respectively. 



domains with sizes of 1, 3, 5 and 7 is lower compared to the model's 
predictions whereas the situation is opposite for domains of size 2, 4, 
6 and 8. Therefore, at equilibrium the DSI model fails to explain the 
saw-tooth shape in the size distribution of charge domains obtained 
using our thermal demagnetization procedure. 

As it is often the case with antiphase boundary formation, a kinetic 
process has to be invoked to explain the observed non-monotonic 
size distribution. Actually, the demagnetization of the system above 
its Curie temperature and its subsequent cooling does not necessarily 
give the system the opportunity to explore the manifold of states and 
to "thermalize". On the other hand, it is established that the system 
does not "remagnetize" homogenously and simultaneously as the 
final configurations reveal long range interactions. Hence, we mod- 
eled the remagnetization process by admitting that "spins" can 
remagnetize completely and sequentially, while interacting with 
already magnetized spins. The interactions are due to both short 
range direct exchange coupling (the elements are connected) and 
long range dipolar interactions. The probability for a spin to remag- 
netize in a given state is directly related to its number of already 
magnetized nearest neighbors and to the energy of the system after 
the possible magnetization (Maxwell- Boltzmann statistics). In other 
words, the remagnetization is driven by the dipolar field generated by 
previously magnetized elements. As represented in figure 3d, the out 
of equilibrium kinetic process we propose yields a behavior quite 
different from the previously presented equilibrium model. Even at 



zero temperature, this finite size system fails to reach the charge 
crystal phase and the charge correlator is limited to —0.69. The 
temperature dependence of the number of domains of small sizes 
also presents some particularities. Although the number of domains 
of size 1, 5, 7 and 8 tends towards zero for small charge correlators, 
the number of domains of size 2, 3, 4 and 6 reachs finite values at 0 K 
and even present a non-monotonic evolution with the charge corre- 
lator (and temperature). Finally, for a given charge correlator, the 
domain size distribution does not decrease regularly, in agreement 
with our experimental findings. The number of domains of size 1 is 
equal to the number of domains of size 2 at a charge correlator of 
— 0.32. The same goes for domains of size 3 and 4 at —0.45 and 5 and 
6 at —0.37. Therefore, a qualitative agreement with our data is found: 
The kinetic model reproduces experimental features that cannot be 
explained by postulating equilibrium conditions. For the sake of 
clarity, the model was kept basic (i.e. just the essential ingredients 
are taken into account to reproduced the main physical observations) 
and we note that the inversion of domain populations occurs at 
charge correlators lower that the one measured experimentally 
( — 0.22). Our main hypothesis was that the magnetization of each 
element switches from "0" to "1" and remains stable. It has been 
recently shown in a system with non- connected elements that the 
elements magnetization can fluctuate at high temperature in a 
"superparamagnetic" manner 101213 . These fluctuations are very effi- 
cient in thermalizing the system. Such fluctuations are not observed 
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at the timescale of our observations but it cannot be excluded that 
they can occur during the cooling down phase. However, the occur- 
rence of such fluctuations is unlikely to give rise to singularities in the 
domain size distribution as it should drive the system towards equi- 
librium. Nevertheless, in order to explore the effect of an incomplete 
thermalization we have run quenched Monte Carlo simulations. A 
limited number of Monte Carlo steps (1-5-25-50-100) are only per- 
formed in order to account for a small number of switching of each 
element. As presented in supplementary information part D, the 
evolution of the charge domains size distribution with such a process 
is qualitatively similar to the equilibrium one. Therefore this model 
of incomplete thermalization fails to reproduce the specific observed 
domain size distribution. 

In conclusion, we find that the remagnetization process of a con- 
nected artificial kagome spin ice leads to a partial order of emergent 
magnetic charges. This observation is compatible with the dipolar 
spin ice model. However a detailed analysis of the charge domain size 
distribution reveals a particular structure (sawtooth distribution). 
This distribution cannot be reproduced by Monte Carlo simulations 
(both at- equilibrium and quenched) and the missing ingredient was 
found to be the kinetics of the remagnetization process. However, it is 
remarkable that a simple out- of- equilibrium process reproduces the 
main features of the equilibrium model. The domain size distribution 
is shown to be a criterion to discriminate between the different pro- 
cesses. These results bring a new insight into the relaxation processes 
of thermally excited artificial spin ice. 
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Methods 

The 2D kagome lattices were fabricated from full films grown by UHV sputtering 
with abase pressure of 10" 9 mbar on Si conducting substrates. Cosputtering of Co and 
Gd in DC mode was used to adjust the concentration of the alloy to get a Curie 
temperature of 475 K. The film composition is Si//Ta(5 nm)/Co 0 7 Gd 0 3 (10 nm)/ 
Ru(2.7 nm). Using e-beam lithography and ion beam etching, we have fabricated 
networks of connected nano-elements with a width of 130 nm. The XMCD-PEEM 
imaging has been performed at the Nanospectroscopy beam line of the Elettra 
Synchrotron (Trieste, Italy) at the L 3 Co absorption edge. 

We performed Monte Carlo simulations for a network of 342 spins that reproduces 
the size and shape of our experimental networks. Long range dipolar interactions 
between spins are taken into account, and the nearest-neighbor coupling is increased 
to account for the finite length of the nanomagnets 8 . The initial temperature was set to 
T/J a p =100 (where J a p is the nearest-neighbor coupling) and corresponds to a high- 
temperature paramagnetic regime. We employed the standard Metropolis algorithm, 
but we also included loop-flips to efficiently drive the system into its ground state. A 
modified Monte Carlo step implies a series of spin and loop flips necessary to achieve 
decorelation between consecutive sampling configurations. The temperature is then 
sequentially reduced to lower energy states until the system reaches the ground state 
manifold. For each temperature value, 10 000 modified Monte Carlo steps are initially 
performed to ensure a good thermalization, followed by 10 000 modified Monte Carlo 
steps used for sampling. The nearest-neighbor charge-charge correlator was com- 
puted at each step and averaged over the entire ensemble for each temperature value. 
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